1 Prepare data

midpoint_rotated <- readRDS("./data/derived/midpoint_rotated.rds")
     midpoint_long <- readRDS("./data/derived/midpoint_long.rds")
midpoint_rotated <- midpoint_rotated %>%
       mutate(
         rate = factor(rate, levels = c("fast", "normal", "slow")),
         stress = factor(stress, levels = c("unstressed", "stressed")),
         item = factor(item),
         V1 = factor(V1, levels = c("a", "e", "i", "o", "u", "y")),
         frame = item,
         duration.log = log(duration),
         z1 = as.numeric(z1),
         z2 = as.numeric(z2),
         f1.z = as.numeric(f1.z),
         f2.z = as.numeric(f2.z)
       )
     
     levels(midpoint_rotated$frame) <- c(
       "b_b", "b_b", "b_w", "b_w", "s_m", "drz_w", "drz_w",
       "k_m", "k_m", "m_m", "m_m", "ni_b", "ni_b", "n_w",
       "n_w", "p_p", "p_p", "p_p", "p_p", "rz_p", "rz_p", 
       "s_m", "z_m", "z_m", "z_w", "z_w"
     )
     
     contrasts(midpoint_rotated$V1) <- "contr.sum"
     
     midpoint_long <- midpoint_long %>%
       mutate(
         rate = factor(rate, levels = c("fast", "normal", "slow")),
         stress = factor(stress, levels = c("unstressed", "stressed"))
       )

2 Normalisation procedure

PL01 <- midpoint_long %>%
       filter(
         speaker == "PL01" &
           fan_line >= 6 &
           fan_line <= 36
       ) %>%
       mutate(
         RateV1 = interaction(rate, V1, sep = "_"),
         StressV1 = interaction(stress, V1, sep = "_"),
         V1 = factor(V1)
       ) 
pl01.gam <- polar_gam(
       Y ~
         V1 +
         s(X) +
         s(X, by = V1),
       data = PL01
     )
## The origin is x = 16.8107470216526, y = -50.1837719761868.
PL01.smooths.plot <- plot_polar_smooths(
       pl01.gam, series = X, comparison = V1
       ) +
       coord_fixed(ratio = 1) + 
       theme_minimal() +
       theme(legend.title = element_blank(),
             text = element_text(size = 6),
             legend.key.size = unit(0.6, "lines")) +
       ggtitle("Tongue contours") +
       scale_x_continuous("tongue position (mm)") +
       scale_y_continuous("tongue height (mm)")
sum_PL01 <- midpoint_rotated %>%
       filter(speaker == "PL01") %>%
       group_by(V1) %>%
       summarize(
         f1.z = mean(f1.z),
         f2.z = mean(f2.z),
         PC1 = mean(PC1),
         PC2 = mean(PC2),
         z1 = mean(z1),
         z2 = mean(z2)
       )
     
     PL01.PCA.plot <- sum_PL01 %>%
       mutate(V1 = factor(V1, levels = c("i", "e", "y", "a", "o", "u"))) %>%
       arrange(V1) %>%
       ggplot() +
         aes(PC1, PC2, label = V1) +
         geom_polygon(color = "black", fill = NA, size = 0.25, linetype = "dashed") +
         geom_text(size = 2.5) +
         scale_x_reverse("PC1") +
         scale_y_reverse("PC2") +
         scale_fill_grey() +
         theme_minimal() +
         coord_fixed() + 
         theme(legend.position = "bottom",
             text = element_text(size = 6),
             legend.title = element_blank()) +
         ggtitle("PCA results")
     
     PL01.rot.plot <- sum_PL01 %>%
       mutate(V1 = factor(V1, levels = c("i", "e", "y", "a", "o", "u"))) %>%
       arrange(V1) %>%
       ggplot() +
         aes(z2, z1, label = V1) +
         geom_polygon(color = "black", fill = NA, size = 0.25, linetype = "dashed") +
         geom_text(size = 2.5) +
         scale_x_reverse("z2") +
         scale_y_reverse("z1") +
         scale_fill_grey() +
         theme_minimal() +
         coord_fixed() + 
         theme(legend.position = "bottom",
               text = element_text(size = 6),
               legend.title = element_blank()) +
         ggtitle("Normalised articulatory space")
     
     PL01.acc.plot <- sum_PL01 %>%
       mutate(V1 = factor(V1, levels = c("i", "e", "a", "o", "u", "y"))) %>%
       arrange(V1) %>%
       ggplot() +
         aes(f2.z, f1.z, label = V1) +
         geom_polygon(color = "black", fill = NA, size = 0.25, linetype = "dashed") +
         geom_text(size = 2.5) +
         scale_x_reverse("F2.z") +
         scale_y_reverse("F1.z") +
         scale_fill_grey() +
         theme_minimal() +
         coord_fixed() + 
         theme(legend.position = "bottom",
               text = element_text(size = 6),
               legend.title = element_blank()) +
         ggtitle("Normalised acoustic space")
ggarrange(
       PL01.smooths.plot, PL01.PCA.plot, PL01.rot.plot, PL01.acc.plot,
       ncol = 2, nrow = 2
     )

ggsave(file = "./figures/norm.pdf", width = 6)
## Saving 6 x 5 in image

3 Descriptive stats plots

3.1 Articulatory undershoot

pl01_rate.gam <- polar_gam(
       Y ~
         RateV1 + 
         s(X) +
         s(X, by = RateV1),
       data = PL01
     )
## The origin is x = 16.8107470216526, y = -50.1837719761868.
plot_polar_smooths(
       pl01_rate.gam,
       series = X,
       comparison = rate,
       facet_terms = V,
       split = list(RateV1 = c("rate", "V")),
       sep = "_"
     ) + 
       coord_fixed(ratio = 1) + theme_minimal() +
       scale_x_continuous("Antero-posterior axis (mm)") +
       scale_y_continuous("Vertical axis (mm)")

pl01_stress.gam <- polar_gam(
       Y ~
         StressV1 +
         s(X) +
         s(X, by = StressV1),
       data = PL01
     )
## The origin is x = 16.8107470216526, y = -50.1837719761868.
plot_polar_smooths(
       pl01_stress.gam,
       series = X,
       comparison = stress,
       facet_terms = V, 
       split = list(StressV1 = c("stress", "V")),
       sep = "_"
     ) + 
       coord_fixed(ratio = 1) + theme_minimal() +
       scale_x_continuous("Antero-posterior axis (mm)") +
       scale_y_continuous("Vertical axis (mm)")

3.2 Duration

means <- aggregate(duration ~  rate, midpoint_rotated, mean)
     means$duration <- round(means$duration, 2)
     
     dur.rate <- ggboxplot(midpoint_rotated, x = "rate", y = "duration") +
       scale_y_continuous("duration (ms)") +
       scale_x_discrete("rate") +
       stat_summary(fun = mean, geom = "point", size = 1) +
       geom_text(data = means, aes(label = duration, y = duration + 5), size = 3) +
       theme(text = element_text(size = 10))
means <- aggregate(duration ~  stress, midpoint_rotated, mean)
     means$duration = round(means$duration, 2)
     
     dur.stress <- ggboxplot(midpoint_rotated, x = "stress", y = "duration") +
       scale_y_continuous("duration (ms)") +
       scale_x_discrete("stress") +
       stat_summary(fun = mean, geom = "point", size = 1) +
       geom_text(data = means, aes(label = duration, y = duration + 5), size = 3) +
       theme(text = element_text(size = 10))
ggarrange(dur.rate, dur.stress, ncol = 2, nrow = 1)

ggsave(file = "./figures/dur_descriptive.pdf", width = 6, height = 3)

3.3 Vowel space

3.3.1 Acoustic vowel space

acc.vowel.plot.rate <- midpoint_rotated %>%
       group_by(V1, rate) %>%
       summarise(mean.f1 = mean(f1.z, na.rm = T),
                 mean.f2 = mean(f2.z, na.rm = T)) %>%
       ungroup %>%
       mutate(V1 = factor(V1, levels = c("i", "e", "a", "o", "u", "y")),
              rate = factor(rate, levels = c("normal", "slow", "fast"))) %>%
       arrange(rate,V1)
## `summarise()` has grouped output by 'V1'. You can override using the `.groups` argument.
acc.rate <- acc.vowel.plot.rate %>%
       ggplot() +
         aes(mean.f2,mean.f1, linetype = rate, label = V1) +
         geom_polygon(color = "black", fill = NA) +
         scale_x_reverse("F2.z", limits = c(2, -1.5)) +
         scale_y_reverse("F1.z", limits = c(2.2, -1.5)) +
         scale_fill_grey() +
         theme_minimal() +
         geom_point() +
         geom_text_repel(size = 4) +
         theme(legend.position = "bottom",
               legend.title = element_blank()) +
         theme(plot.margin = unit(c(0,0, 0,0), "lines"))
     
     acc.vowel.plot.stress <- midpoint_rotated %>%
       group_by(V1, stress) %>%
       summarise(mean.f1 = mean(f1.z, na.rm = T),
                 mean.f2 = mean(f2.z, na.rm = T)) %>%
       ungroup %>%
       mutate(V1 = factor(V1, levels = c("i", "e", "a", "o", "u", "y"))) %>%
       arrange(stress, V1)
## `summarise()` has grouped output by 'V1'. You can override using the `.groups` argument.
acc.stress <- acc.vowel.plot.stress %>%
       ggplot() +
         aes(mean.f2,mean.f1, linetype = stress, label = V1) +
         geom_polygon(color = "black", fill = NA) +
         scale_x_reverse("F2.z", limits = c(2, -1.5)) +
         scale_y_reverse("F1.z", limits = c(2.2, -1.5)) +
         scale_fill_grey() +
         theme_minimal() +
         geom_point() +
         geom_text_repel(size = 4) +
         theme(legend.position = "bottom",
               legend.title = element_blank()) +
         theme(plot.margin = unit(c(0,0, 0,0), "lines"))

3.3.2 Articulatory vowel space

art.vowel.plot.rate <- midpoint_rotated %>%
       group_by(V1, rate) %>%
       summarise(mean.comp1 = mean(z1, na.rm = T),
                 mean.comp2 = mean(z2, na.rm = T)) %>%
       ungroup %>%
       mutate(
         V1 = factor(V1, levels = c("i", "y", "e", "a", "o", "u")),
         rate = factor(rate, levels = c("normal", "slow", "fast"))
       ) %>%
       arrange(rate,V1)
## `summarise()` has grouped output by 'V1'. You can override using the `.groups` argument.
art.rate <- art.vowel.plot.rate %>%
       ggplot() +
         aes(mean.comp2,mean.comp1, linetype = rate, label = V1) +
         geom_polygon(color = "black", fill = NA) +
         scale_x_reverse("z2", limits = c(2, -1.5)) +
         scale_y_reverse("z1", limits = c(1.5, -1.5)) +
         scale_fill_grey() +
         theme_minimal() +
         geom_point() +
         geom_text_repel(size = 4) +
         theme(legend.position = "bottom",
               legend.title = element_blank()) +
         theme(plot.margin = unit(c(0,0, 0,0), "lines"))
     
     art.vowel.plot.stress <- midpoint_rotated %>%
       group_by(V1, stress) %>%
       summarise(mean.comp1 = mean(z1, na.rm = T),
                 mean.comp2 = mean(z2,na.rm = T)) %>%
       ungroup %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "e", "a", "o", "u"))) %>%
       arrange(stress, V1)
## `summarise()` has grouped output by 'V1'. You can override using the `.groups` argument.
art.stress <- art.vowel.plot.stress %>%
       ggplot() +
         aes(mean.comp2, mean.comp1, linetype = stress, label = V1) +
         geom_polygon(color = "black", fill = NA) +
         scale_x_reverse("z2", limits = c(2, -1.5)) +
         scale_y_reverse("z1", limits = c(1.5, -1.5)) +
         scale_fill_grey() +
         theme_minimal() +
         geom_point() +
         geom_text_repel(size = 4) +
         theme(legend.position = "bottom",
               legend.title = element_blank()) +
         theme(plot.margin = unit(c(0,0, 0,0), "lines"))

3.3.3 Save

ggarrange(acc.rate, acc.stress, art.rate, art.stress, ncol = 2, nrow = 2)

ggsave(file = "./figures/descriptive.pdf", height = 6.5, width = 7)

4 Inferential stats

4.1 VIF (multicollinearity)

VIF values below 3 indicate absence of multicollinearity [@zuur2010].

midpoint_rotated %>%
       dplyr::select(duration, f1.z, f2.z, f0.z, rate, stress) %>%
       mutate(
         rate = as.numeric(rate) - 1.5,
         stress = as.numeric(stress) - 1
       ) %>%
       as.data.frame(.) %>%
       usdm::vif(.)
##   Variables      VIF
     ## 1  duration 2.169934
     ## 2      f1.z 1.379401
     ## 3      f2.z 1.208097
     ## 4      f0.z 1.122029
     ## 5      rate 1.639475
     ## 6    stress 1.451903

4.2 Duration

dur.lmer.full <- lmer(
       duration.log ~
         (rate+stress|speaker) +
         (1|frame) +
         rate*V1*stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
ef_dur <- as.data.frame(effect("rate:V1:stress", dur.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("a", "e", "i", "o", "u", "y")))
     
     pd = position_dodge(width = 0.6)
     
     ef_dur <- ef_dur %>%
       mutate(fit = exp(fit),
              lower = exp(lower),
              upper = exp(upper))
     
     dur <- ef_dur %>%
       ggplot() +
         aes(V1, fit, color = stress, linetype = stress) +
         geom_point(position = pd) +
         geom_errorbar(aes(ymin = lower, ymax = upper), position = pd, width = 0, size = 0.5) +
         theme_bw() +
         pl_theme +
         facet_wrap(~rate) +
         scale_x_discrete("Vowel") +
         scale_y_continuous("Duration (ms)") +
         theme(legend.position = "top")
     dur

ggsave(dur, file = "./figures/dur_ef.pdf", width = 7, height = 5)
tab_model(dur.lmer.full)
  duration.log
Predictors Estimates CI p
(Intercept) 3.77 3.63 – 3.91 <0.001
rate [normal] 0.19 0.09 – 0.30 <0.001
rate [slow] 0.45 0.31 – 0.59 <0.001
V11 0.14 0.02 – 0.27 0.023
V12 -0.01 -0.16 – 0.14 0.913
V13 0.00 -0.16 – 0.17 0.968
V14 0.00 -0.12 – 0.13 0.955
V15 -0.22 -0.39 – -0.06 0.008
stress [stressed] 0.23 0.18 – 0.28 <0.001
rate [normal] * V11 -0.01 -0.09 – 0.07 0.786
rate [slow] * V11 -0.04 -0.12 – 0.03 0.268
rate [normal] * V12 -0.01 -0.09 – 0.07 0.780
rate [slow] * V12 0.00 -0.08 – 0.08 0.960
rate [normal] * V13 0.05 -0.03 – 0.13 0.182
rate [slow] * V13 -0.02 -0.10 – 0.06 0.576
rate [normal] * V14 -0.01 -0.08 – 0.07 0.901
rate [slow] * V14 0.02 -0.06 – 0.10 0.593
rate [normal] * V15 -0.04 -0.12 – 0.04 0.308
rate [slow] * V15 0.01 -0.07 – 0.08 0.899
rate [normal] * stress
[stressed] 		0.15 0.10 – 0.20 <0.001
rate [slow] * stress
[stressed] 		0.18 0.13 – 0.23 <0.001
V11 * stress [stressed] 0.12 0.04 – 0.20 0.003
V12 * stress [stressed] -0.05 -0.13 – 0.03 0.215
V13 * stress [stressed] -0.07 -0.15 – 0.00 0.065
V14 * stress [stressed] 0.03 -0.05 – 0.11 0.473
V15 * stress [stressed] -0.01 -0.09 – 0.07 0.841
(rate [normal] * V11) *
stress [stressed] 		-0.03 -0.14 – 0.09 0.655
(rate [slow] * V11) *
stress [stressed] 		-0.06 -0.17 – 0.06 0.331
(rate [normal] * V12) *
stress [stressed] 		-0.04 -0.15 – 0.07 0.465
(rate [slow] * V12) *
stress [stressed] 		-0.04 -0.16 – 0.07 0.455
(rate [normal] * V13) *
stress [stressed] 		-0.07 -0.18 – 0.04 0.209
(rate [slow] * V13) *
stress [stressed] 		-0.08 -0.19 – 0.03 0.173
(rate [normal] * V14) *
stress [stressed] 		0.00 -0.11 – 0.11 0.977
(rate [slow] * V14) *
stress [stressed] 		0.03 -0.08 – 0.15 0.555
(rate [normal] * V15) *
stress [stressed] 		0.12 0.01 – 0.23 0.040
(rate [slow] * V15) *
stress [stressed] 		0.11 0.00 – 0.23 0.046
Random Effects
σ2 0.04
τ00 frame 0.02
τ00 speaker 0.04
τ11 speaker.ratenormal 0.03
τ11 speaker.rateslow 0.05
τ11 speaker.stressstressed 0.00
ρ01 speaker.ratenormal -0.81
ρ01 speaker.rateslow -0.78
ρ01 speaker.stressstressed -0.48
N speaker 10
N frame 12
Observations 1440
Marginal R2 / Conditional R2 0.706 / NA

4.3 F1

f1.lmer.full <- lmer(
       f1.z ~
         (1+stress|speaker) +
         (1|frame) +
         V1*duration*stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
ef1 <- as.data.frame(effect("V1:duration:stress", f1.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "u", "e", "a", "o")))
     
     f1 <- ef1 %>%
       ggplot() +
       aes(duration, fit, linetype = stress, fill = stress) +
       geom_point(data = midpoint_rotated, aes(duration, f1.z, color = stress), alpha = 0.2) +
       geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) +
       geom_line(aes( color = stress)) +
       facet_wrap(~V1) +
       scale_y_reverse("F1.z") +
       scale_x_continuous("duration (ms)") +
       theme_bw() +
       pl_theme +
       theme(legend.position = "top")
     f1

ggsave(f1, file = "./figures/f1_ef.pdf", width = 7, height = 5)
tab_model(f1.lmer.full)
  f1.z
Predictors Estimates CI p
(Intercept) -0.59 -0.78 – -0.41 <0.001
V11 0.38 0.04 – 0.72 0.029
V12 0.08 -0.29 – 0.46 0.660
V13 -0.29 -0.69 – 0.12 0.163
V14 0.31 -0.02 – 0.64 0.064
V15 -0.49 -0.90 – -0.08 0.020
duration 0.01 0.01 – 0.01 <0.001
stress [stressed] 0.40 0.24 – 0.56 <0.001
V11 * duration 0.02 0.01 – 0.02 <0.001
V12 * duration -0.00 -0.00 – 0.00 0.710
V13 * duration -0.01 -0.01 – -0.01 <0.001
V14 * duration -0.00 -0.01 – 0.00 0.180
V15 * duration 0.00 -0.00 – 0.01 0.476
V11 * stress [stressed] 0.41 0.06 – 0.76 0.022
V12 * stress [stressed] -0.22 -0.54 – 0.09 0.165
V13 * stress [stressed] -0.25 -0.55 – 0.05 0.100
V14 * stress [stressed] 0.16 -0.15 – 0.47 0.307
V15 * stress [stressed] 0.18 -0.13 – 0.50 0.246
duration * stress
[stressed] 		-0.00 -0.01 – -0.00 0.005
(V11 * duration) * stress
[stressed] 		-0.01 -0.01 – -0.00 0.024
(V12 * duration) * stress
[stressed] 		0.01 0.00 – 0.01 0.007
(V13 * duration) * stress
[stressed] 		0.00 -0.00 – 0.01 0.608
(V14 * duration) * stress
[stressed] 		0.00 -0.00 – 0.01 0.692
(V15 * duration) * stress
[stressed] 		-0.01 -0.01 – -0.00 0.019
Random Effects
σ2 0.17
τ00 frame 0.07
τ00 speaker 0.00
τ11 speaker.stressstressed 0.01
ρ01 speaker -1.00
N speaker 10
N frame 12
Observations 1440
Marginal R2 / Conditional R2 0.837 / NA

4.4 F2

f2.lmer.full <- lmer(
       f2.z ~
         (1+stress|speaker) +
         (1|frame) +
         V1*duration*stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
ef2 <- as.data.frame(effect("V1:duration:stress", f2.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "u", "e", "a", "o")))
     
     f2 <- ef2 %>%
       ggplot() +
       aes(duration, fit, linetype = stress, fill = stress) +
       geom_point(data = midpoint_rotated, aes(duration, f2.z, color = stress), alpha = 0.2) +
       geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) +
       geom_line(aes( color = stress)) +
       facet_wrap(~V1) +
       scale_y_reverse("F2.z") +
       scale_x_continuous("duration (ms)") +
       theme_bw() +
       pl_theme +
       theme(legend.position = "top")
     f2

ggsave(f2, file = "./figures/f2_ef.pdf", width = 7, height = 5)
tab_model(f2.lmer.full)
  f2.z
Predictors Estimates CI p
(Intercept) -0.06 -0.24 – 0.13 0.550
V11 -0.38 -0.69 – -0.07 0.015
V12 0.28 -0.08 – 0.63 0.132
V13 1.17 0.77 – 1.57 <0.001
V14 -0.52 -0.82 – -0.22 0.001
V15 -0.67 -1.08 – -0.27 0.001
duration 0.00 -0.00 – 0.00 0.283
stress [stressed] 0.05 -0.05 – 0.15 0.332
V11 * duration 0.00 -0.00 – 0.00 0.426
V12 * duration 0.00 -0.00 – 0.01 0.067
V13 * duration 0.01 0.00 – 0.01 <0.001
V14 * duration -0.01 -0.01 – -0.00 <0.001
V15 * duration -0.01 -0.01 – -0.00 0.002
V11 * stress [stressed] 0.04 -0.21 – 0.29 0.765
V12 * stress [stressed] -0.02 -0.25 – 0.20 0.843
V13 * stress [stressed] 0.29 0.08 – 0.51 0.008
V14 * stress [stressed] -0.12 -0.34 – 0.10 0.269
V15 * stress [stressed] -0.26 -0.48 – -0.03 0.025
duration * stress
[stressed] 		-0.00 -0.00 – 0.00 0.600
(V11 * duration) * stress
[stressed] 		-0.00 -0.00 – 0.00 0.630
(V12 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.858
(V13 * duration) * stress
[stressed] 		-0.00 -0.01 – 0.00 0.155
(V14 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.382
(V15 * duration) * stress
[stressed] 		0.00 -0.00 – 0.01 0.566
Random Effects
σ2 0.09
τ00 frame 0.08
τ00 speaker 0.00
τ11 speaker.stressstressed 0.00
ρ01 speaker  
N speaker 10
N frame 12
Observations 1440
Marginal R2 / Conditional R2 0.908 / NA

4.5 Z1

z1.lmer.full <- lmer(
       z1 ~
         (1+stress|speaker) +
         (1|frame) +
         V1*duration*stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
ez1 <- as.data.frame(effect("V1:duration:stress", z1.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "u", "e", "a", "o")))
     
     z1 <- ez1 %>%
       ggplot() +
       aes(duration, fit, linetype = stress, fill = stress) +
       geom_point(data = midpoint_rotated, aes(duration, z1, color = stress), alpha = 0.2) +
       geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) +
       geom_line(aes( color = stress)) +
       facet_wrap(~V1) +
       scale_y_reverse("z1") +
       scale_x_continuous("duration (ms)") +
       theme_bw() +
       pl_theme +
       theme(legend.position = "top")
     z1

ggsave(z1, file = "./figures/z1_ef.pdf", width = 7, height = 5)
tab_model(z1.lmer.full)
  z1
Predictors Estimates CI p
(Intercept) -0.27 -0.41 – -0.13 <0.001
V11 0.17 -0.09 – 0.43 0.211
V12 -0.16 -0.44 – 0.12 0.256
V13 -0.55 -0.85 – -0.26 <0.001
V14 0.25 0.00 – 0.51 0.048
V15 0.12 -0.19 – 0.42 0.453
duration 0.00 0.00 – 0.01 <0.001
stress [stressed] -0.02 -0.15 – 0.11 0.801
V11 * duration 0.01 0.00 – 0.01 <0.001
V12 * duration -0.00 -0.00 – 0.00 0.859
V13 * duration -0.01 -0.01 – -0.00 0.001
V14 * duration 0.01 0.00 – 0.01 0.002
V15 * duration -0.00 -0.01 – 0.00 0.277
V11 * stress [stressed] -0.06 -0.35 – 0.23 0.697
V12 * stress [stressed] 0.07 -0.19 – 0.33 0.591
V13 * stress [stressed] -0.09 -0.34 – 0.16 0.490
V14 * stress [stressed] 0.03 -0.22 – 0.29 0.807
V15 * stress [stressed] 0.04 -0.21 – 0.30 0.737
duration * stress
[stressed] 		-0.00 -0.00 – 0.00 0.619
(V11 * duration) * stress
[stressed] 		-0.00 -0.00 – 0.00 0.722
(V12 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.952
(V13 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.544
(V14 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.891
(V15 * duration) * stress
[stressed] 		-0.00 -0.01 – 0.00 0.591
Random Effects
σ2 0.12
τ00 frame 0.03
τ00 speaker 0.00
τ11 speaker.stressstressed 0.01
ρ01 speaker -1.00
N speaker 10
N frame 12
Observations 1440
Marginal R2 / Conditional R2 0.734 / NA

4.6 Z2

z2.lmer.full <- lmer(
       z2 ~
         (1+stress|speaker) +
         (1|frame) +
         V1*duration*stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
ez2 <- as.data.frame(effect("V1:duration:stress", z2.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "u", "e", "a", "o")))
     
     z2 <- ez2 %>%
       ggplot() +
       aes(duration, fit, linetype = stress, fill = stress) +
       geom_point(data = midpoint_rotated, aes(duration, z2, color = stress), alpha = 0.2) +
       geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) +
       geom_line(aes( color = stress)) +
       facet_wrap(~V1) +
       scale_y_reverse("z2") +
       scale_x_continuous("duration (ms)") +
       theme_bw() +
       pl_theme +
       theme(legend.position = "top")
     z2

ggsave(z2, file = "./figures/z2_ef.pdf", width = 7, height = 5)
tab_model(z2.lmer.full)
  z2
Predictors Estimates CI p
(Intercept) 0.13 -0.11 – 0.37 0.298
V11 -0.27 -0.66 – 0.12 0.182
V12 0.33 -0.14 – 0.79 0.169
V13 1.09 0.57 – 1.62 <0.001
V14 -0.53 -0.92 – -0.15 0.007
V15 -0.59 -1.12 – -0.06 0.030
duration -0.00 -0.00 – -0.00 0.004
stress [stressed] -0.04 -0.15 – 0.06 0.406
V11 * duration -0.00 -0.00 – 0.00 0.650
V12 * duration 0.00 -0.00 – 0.00 0.607
V13 * duration 0.00 0.00 – 0.01 0.002
V14 * duration -0.00 -0.01 – -0.00 0.005
V15 * duration -0.01 -0.01 – -0.00 0.007
V11 * stress [stressed] 0.09 -0.17 – 0.35 0.487
V12 * stress [stressed] -0.09 -0.32 – 0.14 0.466
V13 * stress [stressed] 0.02 -0.20 – 0.24 0.852
V14 * stress [stressed] 0.09 -0.14 – 0.31 0.457
V15 * stress [stressed] -0.19 -0.42 – 0.04 0.104
duration * stress
[stressed] 		0.00 -0.00 – 0.00 0.270
(V11 * duration) * stress
[stressed] 		-0.00 -0.01 – 0.00 0.326
(V12 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.779
(V13 * duration) * stress
[stressed] 		0.00 -0.00 – 0.00 0.658
(V14 * duration) * stress
[stressed] 		-0.00 -0.01 – 0.00 0.312
(V15 * duration) * stress
[stressed] 		0.00 -0.00 – 0.01 0.114
Random Effects
σ2 0.09
τ00 frame 0.15
τ00 speaker 0.00
τ11 speaker.stressstressed 0.00
ρ01 speaker  
N speaker 10
N frame 12
Observations 1440
Marginal R2 / Conditional R2 0.874 / NA

4.7 Composite plot

ggarrange(f2 + theme(legend.position = "top"), z2 + theme(legend.position = "top"), ncol = 2)

4.8 F0

f0.lmer.full <- lmer(
       f0.z ~
         (1+stress|speaker) +
         (1|frame) +
         V1 * duration * stress,
       data = midpoint_rotated,
       control = lmerControl(
         optimizer = "optimx",
         calc.derivs = FALSE,
         optCtrl = list(method = "bobyqa")
       )
     )
     
     summary(f0.lmer.full)
## Linear mixed model fit by REML ['lmerMod']
     ## Formula: f0.z ~ (1 + stress | speaker) + (1 | frame) + V1 * duration *  
     ##     stress
     ##    Data: midpoint_rotated
     ## Control: 
     ## lmerControl(optimizer = "optimx", calc.derivs = FALSE, optCtrl = list(method = "bobyqa"))
     ## 
     ## REML criterion at convergence: 3391.2
     ## 
     ## Scaled residuals: 
     ##     Min      1Q  Median      3Q     Max 
     ## -2.7492 -0.6234 -0.0749  0.5778  4.9360 
     ## 
     ## Random effects:
     ##  Groups   Name           Variance Std.Dev. Corr 
     ##  frame    (Intercept)    0.02512  0.1585        
     ##  speaker  (Intercept)    0.24312  0.4931        
     ##           stressstressed 1.10737  1.0523   -1.00
     ##  Residual                0.54046  0.7352        
     ## Number of obs: 1436, groups:  frame, 12; speaker, 10
     ## 
     ## Fixed effects:
     ##                               Estimate Std. Error t value
     ## (Intercept)                 -0.9271354  0.1902833  -4.872
     ## V11                          0.2990680  0.2423838   1.234
     ## V12                         -0.2139909  0.2340587  -0.914
     ## V13                          0.0721928  0.2234685   0.323
     ## V14                         -0.0418586  0.2258400  -0.185
     ## V15                          0.1287762  0.2420713   0.532
     ## duration                     0.0126725  0.0017110   7.407
     ## stressstressed               0.8131770  0.3577449   2.273
     ## V11:duration                -0.0103272  0.0036658  -2.817
     ## V12:duration                 0.0009472  0.0037718   0.251
     ## V13:duration                -0.0017780  0.0033075  -0.538
     ## V14:duration                -0.0017363  0.0037012  -0.469
     ## V15:duration                 0.0113477  0.0045794   2.478
     ## V11:stressstressed          -0.2810617  0.3197088  -0.879
     ## V12:stressstressed           0.4183847  0.2846774   1.470
     ## V13:stressstressed          -0.2764102  0.2729397  -1.013
     ## V14:stressstressed           0.0117957  0.2794705   0.042
     ## V15:stressstressed           0.3137226  0.2834749   1.107
     ## duration:stressstressed     -0.0085098  0.0020061  -4.242
     ## V11:duration:stressstressed  0.0053035  0.0043042   1.232
     ## V12:duration:stressstressed -0.0075852  0.0044825  -1.692
     ## V13:duration:stressstressed  0.0053459  0.0040823   1.310
     ## V14:duration:stressstressed  0.0003301  0.0042265   0.078
     ## V15:duration:stressstressed -0.0067093  0.0051656  -1.299
## 
     ## Correlation matrix not shown by default, as p = 24 > 12.
     ## Use print(x, correlation=TRUE)  or
     ##     vcov(x)        if you need it
eff0 <- as.data.frame(effect("V1:duration:stress", f0.lmer.full)) %>%
       mutate(V1 = factor(V1, levels = c("i", "y", "u", "e", "a", "o")))
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
     ## predictor f0.z is a one-column matrix that was converted to a vector
eff0 %>%
       ggplot() +
       aes(duration, fit, linetype = stress, fill = stress) +
       geom_point(data = midpoint_rotated, aes(duration, f0.z, color = stress), alpha = 0.5) +
       geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.3) +
       geom_line(aes(color = stress)) +
       facet_wrap(~V1) +
       scale_y_reverse("f0.z") +
       scale_x_continuous("duration (ms)") +
       theme_bw() +
       pl_theme
## Warning: Removed 4 rows containing missing values (geom_point).

ggsave(file = "./figures/f0_ef.pdf", width = 7, height = 5)
## Warning: Removed 4 rows containing missing values (geom_point).
tab_model(f0.lmer.full)
  f0.z
Predictors Estimates CI p
(Intercept) -0.93 -1.30 – -0.55 <0.001
V11 0.30 -0.18 – 0.77 0.217
V12 -0.21 -0.67 – 0.24 0.361
V13 0.07 -0.37 – 0.51 0.747
V14 -0.04 -0.48 – 0.40 0.853
V15 0.13 -0.35 – 0.60 0.595
duration 0.01 0.01 – 0.02 <0.001
stress [stressed] 0.81 0.11 – 1.51 0.023
V11 * duration -0.01 -0.02 – -0.00 0.005
V12 * duration 0.00 -0.01 – 0.01 0.802
V13 * duration -0.00 -0.01 – 0.00 0.591
V14 * duration -0.00 -0.01 – 0.01 0.639
V15 * duration 0.01 0.00 – 0.02 0.013
V11 * stress [stressed] -0.28 -0.91 – 0.35 0.379
V12 * stress [stressed] 0.42 -0.14 – 0.98 0.142
V13 * stress [stressed] -0.28 -0.81 – 0.26 0.311
V14 * stress [stressed] 0.01 -0.54 – 0.56 0.966
V15 * stress [stressed] 0.31 -0.24 – 0.87 0.268
duration * stress
[stressed] 		-0.01 -0.01 – -0.00 <0.001
(V11 * duration) * stress
[stressed] 		0.01 -0.00 – 0.01 0.218
(V12 * duration) * stress
[stressed] 		-0.01 -0.02 – 0.00 0.091
(V13 * duration) * stress
[stressed] 		0.01 -0.00 – 0.01 0.190
(V14 * duration) * stress
[stressed] 		0.00 -0.01 – 0.01 0.938
(V15 * duration) * stress
[stressed] 		-0.01 -0.02 – 0.00 0.194
Random Effects
σ2 0.54
τ00 frame 0.03
τ00 speaker 0.24
τ11 speaker.stressstressed 1.11
ρ01 speaker -1.00
N speaker 10
N frame 12
Observations 1436
Marginal R2 / Conditional R2 0.257 / NA

4.9 F3

midpoint_rotated %>%
       ggplot() +
       aes(duration, f3, linetype = stress, fill = stress) +
       geom_point(alpha = 0.5) +
       geom_smooth(method = "lm") +
       facet_wrap(~V1) +
       theme_bw() +
       pl_theme
## `geom_smooth()` using formula 'y ~ x'